The shape of invasion percolation clusters in random and correlated media 



The shape of two-dimensional invasion percolation clusters are studied numerically for both non- 
trapping (NTIP) and trapping (TIP) invasion percolation processes. Two different anisotropy quantifiers, 
the anisotropy parameter and the asphericity are used for probing the degree of anisotropy of clusters. 
We observe that in spite of the difference in scaling properties of NTIP and TIP, there is no difference in 
the values of anisotropy quantifiers of these processes. Furthermore, we find that in completely random 
media, the invasion percolation clusters are on average slightly less isotropic than standard percolation 
clusters. Introducing isotropic long-range correlations into the media reduces the isotropy of the invasion 
percolation clusters. The effect is more pronounced for the case of persisting long-range correlations. The 
implication of boundary conditions on the shape of clusters is another subject of interest. Compared to 
the case of free boundary conditions, IP clusters of conventional rectangular geometry turn out to be 
more isotropic. Moreover, we see that in conventional rectangular geometry the NTIP clusters are more 
isotropic than the TIP clusters. 

I. Introduction 

Invasion percolation (IP) [TJ [2J [3] is a dynamical percolation process, primarily developed to describe the 
evolution of the interface between two immiscible fluids in a random porous medium. In this process, the 
advance of the interface is modeled as a result of a series of discrete single jumps of the invader (displacing 
fluid) into previously defender (displacing fluid) occupied sites through the least resistant path. The defender 
can be treated as an incompressible fluid. This means that once a portion of it gets surrounded, a trap 
forms and the invader cannot penetrate it further. This variant of invasion percolation is called invasion 
percolation with trapping (TIP). On the other hand, in non-trapping invasion percolation (NTIP) which 
applies for compressible fluids, the invading fluid can potentially enter any region occupied by the defender. 
IP has been also used for modeling corrosion and intrusion 0], simulating the melt infiltration process [5], 
and studying random behaviour of market prices[Hj. In addition to these applications, there are some pure 
scientific interests on the subject. After all, IP is one of the simplest parameter-free models which exhibits 
self-organized criticality [7] [8] . 

Like standard percolation [9], invasion percolation generates self similar fractal clusters. But unlike 
standard percolation, the growth process described above, produces only a single connected cluster. So far, 
much of the efforts have been devoted on investigation of the critical exponents [TOl [HI H2] and scaling 
properties of this cluster [T3] . The statistics of invaded sites and the distribution of sizes of trapped clusters 
in TIP have been studied too [2] [3] Q3] . The shape of IP clusters has remained an open question. 

The shape of random fractals is an important physical property that has been studied for several models 
including lattice animals and percolation clusters [TSl HE1 [T7] , Ising clusters [TS] , random walk [TU] , Eden 
clusters [20] , bond trees [21] and aggregates with tunable fractal dimension [22] . All these studies show that 
anisotropy is an intrinsic property of fractal aggregates. Generally speaking, the shape of a D-dimensional 
cluster is determined by R\ > R\.-- > R 2 D , where Rf's are the eigenvalues (the principal radii of gyration) 
of the cluster radius of gyration tensor 



In the above definition, "af , is the distance of invaded site i from center of mass and N is the size of the 
cluster. If all the R\ are equal, the cluster is spherically symmetric. Otherwise, it is anisotropic and we can 
probe the degree of its anisotropy by defining a proper cluster anisotropy quantifier based on the variations 
in the R 2 [17], which have the following asymptotic form: 
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where v is the leading scaling exponent and is equal to the inverse of D, the fractal dimension of clusters. 
The leading analytic correction-to-scaling term is proportional to iV -1 and N~ e represents the leading 
non-analytic correction-to-scaling term. The coefficients j-j, <Zj, and hi are all independent of TV [15] . 

Two main numerical techniques are commonly used for probing the shape of random clusters. In the 
first method, proposed by Family et al [15] . an asymmetry measure, An = R^/Rf, called the anisotropy 
parameter of an TV-site cluster is evaluated. The quantity An when properly averaged over all clusters with 
the same size is denoted by (An) and is an estimate of the anisotropy parameter of TV-site clusters in the 
ensemble. The case (An) = 1, corresponds to spherical symmetry. For anisotropic objects, (An) is less 
than unity (the term anisotropy parameter may be misleading; the shape of the cluster is more isotropic for 
larger value of An)- The asymptotic behaviour of (A*,) is obtained by taking the limit N — ► oo. Using this 
method for 2-dimension, Family et al, observed for the first time that percolation clusters are not isotropic 
and estimated (A^) = 0.4 as the asymptotic value for the anisotropy of infinitely large percolation clusters. 

The method introduced by Family et al, has this advantage that besides the shape of clusters, it provides 
an un-biased way of evaluating the non-analytical correction-to-scaling exponent [51 US]. Nevertheless, it is 
difficult to treat analytically. A more tractable approach has been suggested by Aronovitz et al [23] and 
Rudnick et al [24] based on the definition of the asphericity Ac as 
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where Q = G R 2 I and R 2 = D^ 1 Tr G. Written in terms of R 2 in 2-dimension, this becomes 
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For an isotropic cluster this quantity is equal to zero. For an ensemble of clusters the asphericity A2 is 
defined to be 
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in which (...) denotes an ensemble average of the quantity. Note that this quantity is different from (A2), 
the ensemble average of A2. Using this method, Quandt et al [18] obtained the value A2 = 0.325 ± 0.006 for 
the asymptotic asphericity of two dimensional percolating clusters, showing again that percolation clusters 
are not isotropic. 

In this paper we study the shape of IP clusters by evaluating both the asphericity and the anisotropy 
parameter. The plan of the work is as follow. After describing the simulation method in section II, we 
present the results of our extensive numerical simulations of the NTIP and TIP processes for completely 
random media in section III. The effect of boundary conditions are examined in section IV. Section V contains 
our estimations of the shape of IP clusters when isotropic long-range correlations are introduced into the 
medium. The paper is concluded at section VI. 



II. Method 



Let us consider a sufficiently large (effectively infinite) square lattice with linear size L, and assign to 
each of lattice sites a random resistance r drawn from an arbitrary distribution D(r). Starting from the 
center of the lattice as a single-site invaded cluster, we follow the growth of the IP cluster by making a series 
of single jumps per time-step to the least resistance neighbor of the cluster. Obviously, the list of the next 
nearest neighbors increases rapidly with time. For the TIP process, we should also consider the possibility 
of formation of traps and discard all the trapped sites from the list of cluster neighbors. In this work, the 
trapping rule has been implemented by using the Hoshen-Kopelmn algorithm |25| . The search for traps is 
time-consuming and makes TIP simulations much slower than NTIP simulations. 

For each cluster of an arbitrary size TV, we evaluate R\ and R 2 , the principal radii of gyration of the 
cluster via diagonalization of the cluster radius of gyration tensor G. The shape of the cluster is then 
characterized by evaluating its asphericity or anisotropy parameter, as described previously. Following the 
growth of the IP cluster in time, we may calculate these values for clusters of any desired size. To achieve 
highly accurate results, we estimate the mean values by sampling the growth of IP cluster in a large number 
of media. The condition of effectively infinite medium requires that none of the IP clusters of a given size 
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N touches any boundary of the medium. More precisely, the linear size of the lattice, L, should be large 
enough, such that all the possible configurations including the most anisotropic ones can potentially appear 
within the lattice boundaries. Otherwise, our sampling will be biased in favor of more isotropic clusters. 

III. The shape of IP clusters in random media 

First we consider the shape of IP clusters in completely random media, i.e. when D(r) is chosen to be a 
uniform distribution. We have followed the growth of IP clusters in 50, 000 different samples and calculate R 2 
and R 2 , f° r a selected values of cluster size in the range 32 < N < 32, 768. The values of N's have chosen such 
that for each block of factor of two in size (e.g. 32 < N < 64, 64 < N < 128,. ..,16384 < N < 32768) there are 
10 equally spaced iV's in the logarithmic scale. For each cluster size N, the anisotropy parameter (An) has 
been calculated by averaging the ratio R\jR\ over different samples. Then, the results have been lumped 
together at the block centers. This procedure not only helps to eliminate correction-to scaling for small 
clusters [18], but it produces new data points which are usually less correlated than the original data [126]. 
The same method has been applied for computing ((R\ — R 2 ) 2 ) and ((iif + R%) 2 ) to obtain the asphericity 
parameter A2 at the centre of each block. The behaviour of anisotropy quantifiers of NTIP clusters are 
depicted in fig.l and fig. 2. For comparison, the anisotropy quantifiers of equilibrium percolation clusters are 
included too. These clusters have been generated using Alexandrowicz method [27] which was later modified 
by Grassberger [55]. In this method, one starts with a single site cluster at the lattice. One of its nearest 
neighbors (perimeter sites) is chosen randomly. This site is occupied with a probability p c — 0.592746, the 
percolation threshold of square lattice. The process continues until the number of perimeter sites becomes 
zero. Only at this point, the radius of gyration tensor is computed. We have generated 400, 000 equilibrium 
percolation clusters of size 32 < N < 32, 768 and compute the ensemble averages (An), ((R 2 — R 2 ) 2 ), and 
((Rf + Rf) 2 ) within each block. 

We observe that when N > 2 10 = 1024 the variation in all curves becomes very small, such that for 
N > 2 12 = 4096, all the curves are effectively flat. This means the effect of correction-to-scaling for 
both NTIP and percolation clusters is negligible and the anisotropy quantifiers have saturated. At this 
limit, the anisotropy parameters of NTIP and percolation clusters fluctuates around 0.337 ± .001, and 
0.389 ± 0.002, respectively. On the other hand, the asymptotic value of the asphericity of NTIP clusters is 
A2 = 0.401±0.002, while for percolation clusters we find A2 = 0.322±0.002. These observations demonstrate 
that NTIP clusters are less isotropic than standard percolation clusters. This is an interesting result, because 
NTIP and standard percolation clusters have the same self-similarity dimension (D — 0.1.8959 ± 0.0001), 
and hence belong to the same universality class [2] [12] . 

How are the An's distributed? To answer this question we have calculated P(A), the normalized distri- 
bution of A for a specified cluster size say, N = 10, 000. To this end, we divided the entire range of [0,1] 
to 50 bins with equal width 5 = 0.02 and counted the number of clusters with the anisotropy parameter 
in the range [^4 — 8, A]. It is seen from fig. 3 that the distribution is asymmetric and quite broad with a 
peak approximately located at A ~ 0.2, which means the most probable configurations are those for them 
R1/R2 ~ V(h2 — 0.45. Our calculation also shows that the fluctuation in A (not shown) is approximately 
equal to 0.19. Furthermore, we observed that the shape of P(A) (and consequently, the fluctuation) is almost 
independent of cluster size N, if N is not too small. 

We have also evaluated the asphericity and the anisotropy parameter of TIP clusters for cluster sizes in the 
range 100 < N < 1000. The results are presented in fig. 4. As it is seen from the figure, there is no difference 
in the shape of TIP and NTIP clusters although the self-similarity dimension of these processes differs from 
each other (D = 1.825 ± 0.005 for TIP in square lattices [T2])- Both (R\) and (R 2 ) have the same leading 
exponents 2v (equation2) and hence, the anisotropy does not involve it. The equivalence of the anisotropy 
quantifiers, therefore, indicates that in addition to the value of a\jai, the ratio of correction-to-scaling terms 
is equal in these processes. 

It is worth to mention that the anisotropy quantifiers are independent of the orientation of the principal 
axes of the cluster, which might be arbitrarily oriented. In fact, the underlying ensembles of clusters are 
isotropic themselves |17) . However, isotropy of an ensemble only implies that a given cluster conformation 
will appear with equal probability in arbitrary orientations |16| . The observed anisotropy in the shape of 
clusters is a result of spontaneous fluctuations in shape about the expected isotropic shape. We may relate 
it to the nature of the dynamics of invasion percolation. As shown by Furuberg et al |3], the advance of 
the interface occurs by invading local areas in bursts; once a new site is invaded, the interface tends to 
stay at that vicinity. Quantitatively, they found that the most probable growth after a time t occurs at a 
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Figure 1: Variation of (An), the anisotropy parameter of equilibrium percolation clusters (diamonds) and 
NTIP clusters (circles) in the range 32 < N < 32, 768. For each factor of two in size, the results have been 
lumped together. The size of errorbars is smaller than the icons used. The absolute value of error in each 
(An) is less than 0.001 for percolation clusters. In the NTIP process, this quantity is less than 0.0005, since 
the number of samples at each block has been much larger than the percolation case. 
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Figure 2: Variation of A2, the asphericity parameter of equilibrium percolation clusters (diamonds) and 
NTIP clusters (circles) in the range 32 < N < 32, 768. The size of errorbars is smaller than the icons 
used. The absolute error in asphericity is not constant and decreases as N increases. It is less than 0.007 
in the NTIP process and slightly greater than 0.001 for equilibrium percolation clusters, when N is large 
(N > 4096). 
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Figure 3: P(A), the normalized distribution of A for NTIP clusters of size N = 10, 000. 

distance d t ~ i 1 / 2 , where z is the dynamic exponent. Naturally, this local growth might amplify any small 
fluctuations in the ratio of R\jR\. 

IV. The shape of IP clusters in conventional geometry 

In the more conventional simulations of invasion percolation processes, the host lattice is assumed to be 
a L x 2L rectangular lattice, and instead of the center, the invasion process starts from one of the smaller 
lattice edges. The outlet or sink is located on the opposite edge and the other two lattice edges are assumed 
to be impermeable. The growth process stops at breakthrough, when the invader reaches the outlet. In this 
situation, the IP cluster connects the inlet and outlet through a single, continuous path. The properties 
of this sample spanning cluster(SSC) within the central L x L part of the lattice, i.e. far from inlet and 
outlet [3J, have been the subject of intense research. 

To estimate the asymptotic value of the anisotropy parameter of the central part of SSC, we generated 
10, 000 samples for each of lattice sizes, L = 64, L = 128, L = 256, and 2000 samples of size L = 512. The 
mean anisotropy parameter Al is then computed for each L. In this geometry, the mass of SSC varies in 
different realizations even when L is fixed. For example, in ordinary TIP the mass of central part of SSC is 
N = (5.4 ± 1.1) x 10 4 for L = 512. Nevertheless, since iV is very large itself, this variation does not affect 
the value of A^ via correction-to-scaling terms. In fact, our simulations show that Al does not depend 
on L, if L is sufficiently large. The obtained value of Aoo is 0.64 for the NTIP process, and 0.57 for TIP 
process. Compared to the previous case, the shape of SSC in both NTIP and TIP has turned out to be 
more isotropic. This is because in this case, the growth process continues even after the IP cluster touches 
the boundaries of the central L x L frames. The difference between the shape of A^ in this geometry is a 
consequence of trapping rule which limits the growth of SSC in the TIP process. 

V. The effect of long range correlations on the shape of IP clusters 

In many practical applications, the nature of disorder is not completely random and there are correlations 
in the properties of the medium [521 [3D]. To investigate the effect of correlations on the shape of IP clusters, 
we have considered the case for which the distribution of the resistance of lattice sites obeys the statistics 
of fractional Brownian motion (FBM) Z?s(x) [3TJ [34] . FBM is a stochastic process whose increments are 
statistically self-similar such that its mean square fluctuation is proportional to an arbitrary power of the 
spatial displacement x 
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Figure 4: Comparison between the asphericity (upper curve) and the anisotropy parameter(lower curve) of 
TlP(circles) and NTIP (dots) clusters. 



([D B (x)-D B (0)f)~\ X \ 2H (6) 

H is called the Hurst exponent and determines the type of correlations. If H = 0.5, the above equation pro- 
duces the ordinary Brownian motion, which means that in this case there is no correlation between different 
increments. If H > 0.5, then FBM generates positive correlations, i.e. all the points in a neighborhood of a 
given point obey more or less the same trend. If H < 0.5, FBM is anti-persistence, i.e. a trend at a point 
will not be likely followed in its immediate neighborhood. 

The reason that we have chosen FBM process is twofold. First, FBM generates long-range and at the same 
time isotropic correlations in the field. Therefore, the host lattice retains its isotropy. Second it has been 
demonstrated that such process has practical applications in earth sciences and also reservoir engineering, 
where the permeability field and also the porosity distribution of many real oil reservoirs and aquifer follow 
FBM statistic (S3 E21 E3] . 

There are a number of methods which are capable of producing the FBM statistics j3H GH] ■ We have 
used one of the most popular one, the method of fast Fourier transformation (FFT) filtering which is based 
on the fact that the power spectrum of FBM is given by: 

where is a numerical constant, to = (u>x, lu 2 ), with uii being the Fourier component in the ith direction and 
p = H + 1. In FFT method, one starts with a white noise W(x,y) defined on the lattice sites. The power 
spectrum of W(x, y) is constant and independent of frequency. Therefore, filtering W(x,y) with a transfer 
function yj Sb{u) generates another noise whose spectral density is proportional to Sb(lu). The method is 
straightforward and fast, but it usually produces periodic noises. Therefore, one has to produces a larger 
lattice and keeps only a portion(typically j in two dimensional lattices). 

In fig. 5 we have reported our estimation of the anisotropy parameters of IP clusters in media obeying 
the FBM statistics in the range 10 < N < 800. In this figure, we have compared the value of An for 
three different Hurst exponents, H = 0.2 (anti-persistent correlation), H = 0.8 (persistent correlations), and 
H = 0.5 (Brownian motion) with the results of completely random media. The data have been obtained 
from averaging over 30, 000 samples for each case. Like completely random media, we observed no difference 
between the shape of NTIP and TIP clusters (not shown). These results indicate that any deviation from 
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Figure 5: The anisotropy parameters of IP clusters in media obeying the FBM statistics; H — 0.2 (diamonds); 
H = 0.5 (squares), and H = 0.8 (triangles). For comparison the data for completely random media (circles) 
has been included too. 



complete randomness makes the shape of invasion percolation clusters more anisotropic. Furthermore, we 
find that IP clusters in the presence of persistent correlations are less isotropic than IP clusters of anti- 
persistent correlations. Based on what has been explained in the last lines of section III, these effects can be 
assigned to the difference between dynamics of invasion percolation in random and correlated media. In fact, 
we anticipate that the burst-like growth occurs more effectively, maybe with different dynamic exponent and 
amplitude (which depend on the nature of the disorder), resulting more anisotropy in the shape of clusters. 
The difference between the shape of clusters for H = 0.2 and H = 0.8 is compatible with this image. The 
presence of persistent long-range correlations intensify the burst-like growth and as the result, IP clusters 
become more anisotropic in this case. 

VI. Conclusions 

The shape of clusters in IP processes have been probed numerically by evaluating their asphericity and 
anisotropy parameters. The results indicate that the shape of clusters are the same for both TIP and NTIP 
processes. This conclusion does not depend on the type of disorder in the host lattice. We found that similar 
to other random fractals, generated in a variety of stochastic processes, the invasion percolation clusters are 
anisotropic too. Moreover, we observed that IP clusters are less isotropic than standard percolation clusters. 
By introducing long-range correlation into the media the clusters became more anisotropic in shape than 
before. These effects might be explained according to the dynamics of invasion percolation and the burst-like 
nature of the growth process of IP clusters. 
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